Viscoelasticity near the gel-point: a molecular dynamics study 
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We report on extensive molecular dynamics simulations on 
systems of soft spheres of functionality /, i.e. particles that 
are capable of bonding irreversibly with a maximum of / other 
particles. These bonds are randomly distributed throughout 
the system and imposed with probability p. At a critical 
concentration of bonds, pc ~ 0.2488 for / = 6, a gel is formed 
and the shear viscosity 77 diverges according to ~ (pc— p)""- 
We find s ~ 0.7 in agreement with some experiments and with 
a recent theoretical prediction based on Rouse dynamics of 
phantom chains. The diffusion constant decreases as the gel 
point is approached but does not display a well-defined power 
law. 
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The behavior of transport coefficients and elastic mod- 
uli near the gelation transition has been discussed in the 
literature for many years [Q. To date no consensus on 
either the theoretical or experimental side has emerged 
as far as the critical behavior of these quantities is con- 
cerned. The phenomenology is as follows. As monomers 
or polymers are randomly crosslinked to each other in a 
melt, the shear viscosity 77 increases with crosslink con- 
centration p and diverges at a critical concentration pc 
at which an amorphous rigid network is formed. Exper- 
iment and theory both yield 77 ^ {Pc ~ p)~'^ but there 
is no general agreement regarding the value of the expo- 
nent s. Indeed there is good reason to expect at least two 
different universality classes: As de Gennes showed 
vulcanization (crosslinking of very long chains) must be 
distinguished from gelation (crosslinking of short chains 
or monomers) as far as critical behavior is concerned. 
This conclusion is supported by recent experiments 
which show quite clearly that chain length is a relevant 
parameter. 

Before describing our model and calculations, we dis- 
cuss briefly the experimental and theoretical situation 
for gels putatively in the percolation (short chains) uni- 
versality class. One group of experiments has produced 
exponent values for the shear viscosity in the range 
0.6 < s < 0.9 §,§. Another group has reported 

values for s in the range 1.1 — 1.3 and interpreted ^] 
these in terms of a Rouse model without hydrodynamics. 
On the theory side, de Gennes |^ argued that the viscos- 
ity is analogous to the conductivity of a random mixture 
of normal conductors and superconductors, with an ex- 



ponent s w 0.67. The aforementioned Rouse model g] 
prediction is s = — /3 « 1.35, where v « 0.88 and 
f3 « 0.41 are the correlation length and order parame- 
ter exponents of percolation theory in three dimensions. 
Finally, a recent theoretical analysis of a different model 
with Rouse-dynamics [|l^ has predicted s — 4> — (3 ~ 0.7 
where (j) ~ l-H is the crossover exponent of a random 
resistor network. Given this wide disparity in both the- 
oretical and experimental results, computer simulations 
may help to clarify the situation. 

The shear modulus /i of a rigid network near the gel 
point is typically entropic in nature and it vanishes with 
a power law of its own as the gel point is approached from 
the rigid phase /i ^ {p — PcY . Recent numerical work on 
systems in the percolation universality class has pro- 
vided evidence that t — f in both two and three dimen- 
sions, where / is the exponent that describes the critical 
behavior of a randomly disordered network of conductors 
and insulators near the percolation point. This result is 
consistent with another argument of de Gennes p^ . In 
the dynamical scaling theory of the gelation transition, 
the two exponents t and s are not independent, but rather 
obey the sum rule s + 1 ^ z where z describes the diver- 
gence of the longest relaxation time in the incipient gel: 
t* — tQ{pc — p)^^ |l). This connection allows an impor- 
tant consistency check between the results reported here 
and those of [ pT[ . 

The model that we simulate is capable of describing 
the entire range from simple liquid to entropic solid. 
All particles interact through a soft sphere potential 
ysa{rij) — e{a/rijY^ [|l3| with a — 1 and, for our sim- 
ulations, kBT/e = 1. If there are no other interactions, 
this system forms a simple three-dimensional liquid at 
least at low density. All of our simulations are done at a 
volume fraction $ = iiNa^ /QV = 0.4 which is well away 
from the liquid-solid coexistence density. The viscosity of 
the system is progressively increased by introducing ran- 
dom crosslinks between particles. Specifically, the system 
of particles is initially placed on a simple cubic lattice 
that fills the cubic computational box. Each particle may 
bond with probability p with each of its six nearest neigh- 
bors. This bonding is permanent and enforced by the 
spherically symmetric potential Ynnifij) = \k{rij — tq)^ 
with k — Se/cr^ and ro = (7r/6<I>)^/'^cr. For p < pc, where 
Pc ~ 0.2488 is the bond percolation probability of the 
simple cubic lattice, the system consists of finite clusters 
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of varying masses. For p > pc the system is an entropic 
solid with nonvanishing shear modulus. This system has 
been previously studied by us for Vsa = and p > pc- 
Farago and Kantor and Cohen and Plischke [Q have 
shown that self-avoidance is irrelevant as far as the crit- 
ical behavior of the elastic constants is concerned. 

In the simulations, the system of particles is first equi- 
librated for 5x10^ time steps with Brownian dynamics. 
At the end of this equilibration time, the damping and 
thermal noise are turned off and the subsequent evolution 
is conservative. The equations of motion are integrated 



with a standard velocity Verlet algorithm 1 15 1 with a time 
step of St = O.OOS-y/mcr^/e. The shear viscosity 7y(p) and 
the self-diffusion constant D{p) are then calculated from 
the appropriate Green-Kubo formula |lq]: 



rj = lim / dtCcra (t) 
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and where V is the derivative of the total potential en- 
ergy. Similarly, the diffusion constant is given by the 
familiar expression: 
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It is well known that the velocity-velocity correlation 
function decays extremely slowly, typically with a 'long 
time tail' t~^^'^ power law, even in simple liquids 
We find the same long-time behavior in our simulations 
as well, remarkably for all values of p. There is consid- 
erably more disagreement regarding the stress-stress cor- 
relation function. Extended mode-coupling theory jT^ ] 
suggests that close to the melting point the stress cor- 
relator decays exponentially at long times. Powles and 
Heyes have found that both an exponential decay 
and a Lorentzian provide a reasonably good fit to their 
data in the dense liquid regime. In our situation where 
clusters of various sizes form the system, the decay is 
dramatically affected by crosslinking and becomes very 
slow close to the gel-point. Therefore in the evaluation 
of (|l|,||) we have used time series from tmax = 300r to 
tmax = 1200r where r is the average time between col- 
lisions for particles in the uncrosslinked liquid (j> = 0). 
Even with such long runs, for p close to Pc an estimate 
of the residual integral to i = cxd had to be added. This 
is discussed further below. 
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FIG. 1. The dimensionless stress-stress correlation function 
(jCaa/rn for L = 12 and several crosslinking probabilities for 
(a) all t/r > and (b) for t/r > 2. In the uppermost curve 
in (b) we also show the fit of Caa to a stretched exponential 
(solid line). For p — 0.1 the data are obtained from 30 dif- 
ferent crosslinkings, for larger p from between 100 and 200 
crosslinkings. 

We have simulated systems consisting of = L'^ par- 
ticles with L — 5,8, 12 and 20, the first three over the 
concentration range < p < 0.24, the fourth only for 
p > 0.20. For these relatively small systems, the proba- 
bility that one of the clusters percolates in at least one 
of the three directions is an issue. If there is perco- 
lation, then the viscosity is not a well-defined quantity 
and the sample is characterized instead by a shear mod- 
ulus. Therefore, we have eliminated all percolating sam- 
ples from the calculation. For L = 5, a non negligible 
number of samples percolates in at least one direction 
a,t p — 0.15; at L = 20 percolation becomes significant 
only at p — 0.23. Figure |l| depicts the stress correlator 
for iV = 12^^ particles for a set of crosslink concentra- 
tions. The top panel shows this function over the entire 
range in time, the lower panel for t > 2r. The top panel 
shows that the changes for short times between the liq- 
uid {p = 0) and the incipient gel {p = 0.24) are extremely 
small. The effect of increasing crosslink density is illus- 
trated in the lower panel where it is clearly seen that the 
decay of the correlation function becomes progressively 
slower as the gel point is approached and that even at the 
longest time, for p > 0.2 the correlator is nonnegligible. 
Therefore, integrating Caa(t) only to tmax would result in 
an underestimate of the shear viscosity. In order to cap- 
ture the remaining contribution, we have fit Caa (t) with 



2 



a stretched exponential over various windows {ts,tmax) 
for starting values tg > 2t. Such a fit is shown in Figure 
|T|b for the uppermost curve {p = 0.24) over the range 
2t < t < trnax- The fit is essentially indistinguishable 
from the simulation. However, there is no fundamen- 
tal reason to believe that a stretched exponential is the 
functional form of the long-time behavior of Co-a- (t) . 
This is the reason that we have chosen several different 
starting points for the fit: the spread in values of the 
remaining integral of the fitting function from tmax to 
infinity provides an estimate of the error associated with 
this part of the calculation. 

In Figure |2| we display our data for the dimensionless 
shear modulus r/a^ /^/mkBT in two versions. In part (a) 
the raw data is shown as function of {pc — p) together 
with a guide to the eye of the form a{pc ~ p)"'^. This 
function clearly captures the general behavior of the data 
in the intermediate range of p. For p close to zero, one 
would not expect the system to anticipate the formation 
of a gel at p « 0.25 and for p close to Pc finite-size ef- 
fects are clearly evident. Part (b) of this figure attempts 
to collapse the data by means of the finite-size scaling 
ansatz ??(L,p) = L^/'''^[L(j) - pc)"] with s = 0.7 and 
V = 0.876. Internal consistency requires that the scal- 
ing function have the asymptotic form ^'(a;) ^ x~''l^ for 
large x and a line corresponding to this form is also shown 
on the figure. 
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Further support for the conclusion s « 0.7 comes 
from the complex frequency-dependent viscosity "rfii^^ = 
r({uj) +iri"{u)) = G*{Lu)/{iuj) where G* is the complex 
elastic modulus. The scaling ansatz for these functions 0| 
is lim^^_^o v'i^) ~ iPc - p)^" for p < Pc, lim^^o G'{uj) ~ 
{p — PcY for p > Pc and, for frequencies co > uj* , where 
is a characteristic crossover frequency that approaches 
zero as p ^ pc, G*{lu) ~ (iuj)". The connection between 
the critical behavior of the modulus in the rigid phase and 
the viscosity in the fluid phase is then expressed through 
the scaling relation u = t/{s + 1). Moreover, in the high 
frequency region, one expects r]'{uj) and rj"{ijj) to both 
vary as o;"^^ and the ratio of the real and imaginary parts 
to obey u — 2/Trtaii~^{ri' /ij"}. In our previous work on 
the rigid phase , wc have concluded that t « 2 in three 
dimensions. Therefore, with s « 0.7 we have u « 0.74. 
The frequency-dependent viscosity is plotted in Fig. ^ 
for L = 12 and p > 0.20. There is clearly a region of 
power-law behavior that extends to lower frequencies as 
the critical point is approached. This behavior is seen 
more clearly in rj' than in 77". Nevertheless, both pieces 
of the shear viscosity decrease in a way consistent with 
j^~o.27 jj-^ very satisfactory agreement with the foregoing 
analysis. As well, the ratio of rj' to rj" in this regime 
produces a second estimate u w 0.76. 




FIG. 3. Plot of the complex viscosity r]*{Lj) as function of 
LO for L = 12 and p = 0.2 (lowest curve in each set), 0.22, 
0.23 and 0.24 (top curve in each set). The power law form 
rj ~ oj""^ is more evident for rj' than for 77" and becomes more 
prominent as p — > pc. 
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FIG. 2. Log-log plot of the dimensionless shear viscosity 
as function of crosslinking probability p: (a) raw data; (b) 
finite-size scaling form of the data. 



Finally, in Figure 
diffusion constant 



we show the dimensionless self- 
obtained 



m/ (j'^kBTD obtained from the 
velocity- velocity correlator (||) for all particles in the sys- 
tem. In this case, we show only the raw data. It seems 
clear from the behavior of D in the critical region that 
a finite-size scaling analysis is unlikely to improve the 
collapse of the data. For p < 0.2 the data are not incon- 
sistent with a power law behavior of the form {pc — p)°'^ 
but the evidence for this is weak at best. Moreover, the 
fact that the data for p ~ Pc are essentially independent 
of L suggests that D{L — > 00, p Pc) is finite. Pre- 
cisely at the gel-point, in the thermodynamic limit, the 
system consists of a percolating cluster with fractal di- 
mension Dp « 2.5. The particles that are not on the 
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infinite cluster are organized into finite clusters of vari- 
ous sizes. Approximately 18% remain as monomers that 
presumably are able to diffuse quite easily through the 
percolating cluster since this cluster contains holes on 
all length scales. This would account for the absence of 
critical behavior in D. 
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FIG. 4. The dimensionless diffusion constant D{L,p) as 
function of p. The straight line corresponds to a power law 

We are aware of one previous simulation that at- 
tempted to address the critical behavior of the shear vis- 
cosity near the gel point. Recently del Gado et al. | |T9| ] 
simulated a very different model, namely particles con- 
fined to the sites of a lattice and randomly crosslinked 
to form clusters of various sizes. This system was then 
evolved by a bond fluctuation method and the diffusion 
constants Dm of clusters of mass m measured. They pos- 
tulated the relation D{R) - R-^^+''/''\ where R is the 
radius of gyration of a cluster, and from this determined 
s ~ 1.3. We are not aware of any rigorous derivation 
of this connection between diffusion and viscosity. How- 
ever, it may be that their model simply contains different 
physics. 

In conclusion, we have obtained the shear viscosity of a 
viscoelastic liquid below the gel point for a conceptually 
simple model that we have previously studied in the rigid 
phase. We have presented evidence that the shear vis- 
cosity diverges with an exponent s « 0.7 consistent with 
recent theoretical work |]l0| and some of the experimental 
data 1^. Future extensions of this work will include an 
investigation of viscoelasticity in two dimensions and a 
study of dynamics in the gel phase. 
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